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Abstract 

Reaction networks are systems in which the populations 
of a finite number of species evolve through predefined 
interactions. Such networks are found as modeling tools 
in many disciplines, spanning biochemistry, epidemiology, 
pharmacology, ecology and social networks. It is now 
well-established that, for small population sizes, stochastic 
models for reaction networks are necessary to capture ran- 
domness in the interactions. The tools for analyzing them, 
however, still lag far behind their deterministic counter- 
parts. In this paper, we bridge this gap by developing 
a constructive framework for examining the long-term be- 
havior and stability properties of the reaction dynamics in 
a stochastic setting. In particular, we address the prob- 
lems of determining ergodicity of the reaction dynamics, 
which is analogous to having a globally attracting fixed 
point for deterministic dynamics, and determining mo- 
ment bounds for the underlying stochastic process. Theo- 
retical and computational solutions for these problems are 
obtained by utilizing a blend of ideas and techniques from 
probability theory, linear algebra, polynomial analysis and 
optimization theory. We demonstrate that stability prop- 
erties of a wide class of networks can be assessed from the- 
oretical results that can be recast as efficient and scalable 
linear programs, well-known for their tractability. It is 
notably shown that the computational complexity is often 
linear in the number of species, but worst-case quadratic. 
We illustrate the validity, the efficiency and the universal- 
ity of our results on several reaction networks arising in 
fields such as biochemistry, epidemiology and ecology. 



1 Introduction 

Reaction networks represent a modeling paradigm that 
finds applications in many areas of science. Examples in- 
clude, chemical reaction networks [S], cell signalling net- 
works [28], gene expression networks [32], metabolic net- 
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works [30] J pharmacological networks [3], epidemiologi- 
cal networks [T7j and ecological networks [2]. Tradition- 
ally, reaction networks are mathematically analyzed by 
expressing the dynamics as a set of ordinary differential 
equations. Such a deterministic model is reasonably ac- 
curate when the number of network participants is large. 
However, when this is not the case, the discrete nature of 
the interactions becomes important and the dynamics is 
inherently noisy. This random component of the dynam- 
ics cannot be ignored as it can have a significant impact 
on the macroscopic properties of the system [121 IMl [23] • 
To account for this randomness and study its effects, a 
stochastic formulation of the dynamics is necessary. The 
most common approach is to model the dynamics as a 
continuous-time Markov process whose states denote the 
current population size. Many recent articles have been 
devoted to stochastic models and their analysis in view 
of understanding the role of noise and its impact on the 
system's behavior [8l IT] [22l [25] . 

Among the two modeling approaches, the determinis- 
tic approach is far more tractable than the stochastic ap- 
proach. This is not surprising since our knowledge of or- 
dinary differential equations is fairly advanced. For deter- 
ministic models, many tools allowing us to analyze the sys- 
tem without simulating the trajectories are indeed avail- 
able. Unfortunately, this is not the case for stochastic 
models. Most papers that use such models rely on heavy 
computational techniques such as the simulation of several 
trajectories in order to determine the relevant character- 
istics of the system. More importantly, since one can only 
simulate a finite number of trajectories for a finite amount 
of time, properties like long-term behavior and stability 
cannot be verified through such simulations. In this pa- 
per, we overcome this problem and find a direct way to 
examine such properties for the stochastic model, without 
relying on simulations. Our approach combines stochas- 
tic analysis along with linear algebra, polynomial analysis 
and optimization techniques. For the latter we restrict 
ourselves to Linear Programming in this paper. Refine- 
ments of our method using more advanced optimization 
techniques will be considered elsewhere. 

2 Reaction networks 

Before we introduce the properties we study in this pa- 
per, we formally describe our reaction network. Moti- 



vated by the literature on chemical kinetics, we refer to 
the network participants as molecules which may belong 
to one of d species. There are K reactions in the net- 
work and for any k = I, . . . , K, the stoichiometric vector 
Ck = (Cfc.ii ■ ■ • 1 Ck,d) denotes the change in the number of 
molecules in each of the species due to the fc-th reaction. 
When the state of the system is x, the fc-th reaction fires 
at rate Xk{x). The functions Ai, . . . , Xk are known as the 
propensity functions in the literature. The notion of state 
is different for deterministic and stochastic approaches. In 
the deterministic setting, the state is a vector of concen- 
trations of the d species, while in the stochastic setting, 
the state refers to the vector of molecular counts of the d 
species. 

2.1 Deterministic models 

Consider the deterministic model for the reaction network 
described above. If the initial state is xq, then the evolu- 
tion of concentrations is given by {4'xoit))t>o which satis- 
fies the Reaction Rate Equations (RRE) of the form 



d(t)xoit) 
dt 



K 

E 
fe=i 



Afe(0xo(*))Cfe with (j)xo{Q)^XQ. (1) 



We are interested in the long-term behavior and stability 
of our reaction dynamics. More precisely, we would like 
to check if the following conditions are satisfied. 

DCl For any xt^, there is a compact set /C(xo) such that 
(j)xa{t) elC{xo) for ah t > 0. 

DC2 There exists a compact set /Cq such that for any xq, 
we have (f>xg {t) G /Co for large values of t. 

DCS There is a x^q such that for any xq we have (Jj^q (t) — )■ 
Xcq as t — >■ oo. 

The first condition, DCl, says that for any Xq, the en- 
tire trajectory {(l)xoit))t>o stays within some compact set. 
We would expect this to be true for most realistic sys- 
tems. Hence a violation of this property may suggest a 
flaw in the deterministic model. The second condition, 
DC2, says that there is an attractor set for the dynam- 
ics, where all the trajectories eventually lie, irrespective 
of their starting point. This is related to the notion of 
permanence for reaction networks. The last condition, 
DCS, says that there is a globally attracting fixed point 
for the deterministic model. Given a reaction network one 
can verify these conditions by using techniques from the 
theory of dynamical systems (211 I31j , without having to 
simulate the paths of the deterministic process. There is 
also a general theory available that aims to check con- 
dition DCS for reaction networks satisfying mass- action 
kinetics (see [TTl [TH ^l HHj). Our goal in this paper is 
to develop a theoretical and computational framework for 
verifying conditions similar to DC1,DC2 and DCS for 
stochastic models of reaction networks. 

2.2 Stochastic models 

Consider the stochastic model corresponding to the reac- 
tion network described above. In this setting the firing 
of reactions are discrete events. When the state of the 



system is x, the fc-th reaction fires after a time which is 
an exponentially distributed random variable with rate 
Xkix). The dynamics can be represented by the Markov 
process {X^g (i))t>o where Xq is the initial state. Note that 
a XxfXt) — (Xi{t), . . . ,Xd{t)), then Xi{t) is the number 
of molecules of the i-th species at time t. It is important to 
choose a suitable state space S for this Markov process. If 
No is the space of non-negative integers then let S be the 
smallest non-empty subset of Ng satisfying the following 
property: if x G 5 and Xk{x) > for some k = 1, . . . , K, 
then x + (k G S. Observe that ii xq E S then we also have 
that Xxo{t) e S for all i > 0. Hence S can be taken to 
be the state space of all the Markov processes described 
as above with the initial state in S. 

Let Vis) denote the space of probability distributions 
over S, endowed with the weak topology (see PH]). A set 
/C C V{S) is compact in this topology if for any e > 
there is a finite set A^ C S such that 

sup ^(Ag) > 1 — e. 

For any x,y E S let px{t,y) denote the following proba- 
bility 



Px{t,y)^V{Xx{t)^y). 



Note that 



Px(0,y) 



1 \i X ^ y 
otherwise. 



(2) 
(3) 



Defining Px{t,A) = Y.yaAPx^^-'V) ^^^ ^^y A C 5, we can 
view Px{t) as an element in V{S). In fact, Px{t) is the dis- 
tribution at time t of the Markov process {Xx{t))t>o- The 
dynamics of Px (t) is given by Kolmogorov's forward equa- 
tion which is also called the Chemical Master Equation in 
the literature. It has the following form: 



K 



K 



-^^ — = ^Pxit,y - Ck)Xkiy - Ck) -Px{t,y)^Xk{y). 

fe = l k = l 

(4) 



Theoretically, one can find Px {t, y) for any i > and y E S, 
by solving this system with initial condition ([3]) . However 
this system consists of as many ordinary differential equa- 
tions as the number of elements in S. This system can 
be explicitly solved if S is finite, but this only happens in 
very restrictive cases where all the reactions preserve some 
conservation relation. Typically, S is infinite and solving 
this system analytically or even numerically is quite dif- 
ficult, if not impossible. Our results are not very useful 
when S is finite. For most problems of practical interest 
S is infinite and this is what we assume from now on. 

The above discussion shows that at the level of distri- 
butions, we can view the stochastic dynamics {Xxg{t))t>o 
as the deterministic dynanfics {Pxo{'t))t>o^ which satisfies 
the system Q. However, the major difficulty in analyz- 
ing this deterministic dynamics is that it occurs over an 
infinite dimensional space V{S). Nevertheless we can re- 
cast the conditions DCl, DC2 and DCS in the stochastic 
setting as below. 

SCI For any xq, there is a compact set JC{xo) C 1^(5) 
such that Pxo{t) S /C(a;o) for all i > 0. 



SC2 There exists a compact set /Co C 7^(5) such that for 
any xq G 5 we have Pxoit) G ^o for large values of t. 

SC3 There is a tt e VIS) such that for any xq we have 

PxQ (t) — > TT as t — > oo. 

Each of the above conditions give an important insight 
about the long-term behavior and stability of the stochas- 
tic dynamics. The first condition, SCI, says that for every 
e > we can find a finite set A^ C S such that each p^g (t) 
puts at least (1 — e) of its mass in A^. In other words, 
for any t > the probability that X^o^t) hes outside A^ 
is less than e. This also means that the Markov process 
representing the reaction network remains bounded with 
time, which should be true for most realistic models. If 
condition SC2 holds then the evolution of distributions 
have a compact attractor set in 'P{S), where all the tra- 
jectories eventually lie irrespective of their starting point. 
This suggests that in the long run, the family of processes 
{{Xxo(t))t>o '■ xq G S}, spend most of their time on the 
same set of states. The last condition SC3 says that the 
evolution of distributions have a globally attracting fixed 
point TT. If this holds, then the Markov process represent- 
ing the reaction dynamics is ergodic with n as the unique 
stationary distribution. When it comes to understanding 
the long-term behavior of a stochastic process, then er- 
godicity is a desirable property to have. In the long-run, 
the proportion of time spent by any trajectory of an er- 
godic process, in any subset of the state space is equal 



form 



to the stationary probability of that subset (see (111). In 
other words, information about the stationary distribu- 
tion can be obtained by observing just one trajectory for 
a sufficiently long time. Such a result is very important 
in applications. For example, consider a culture with a 
large number of identical cells with each cell having the 
same reaction network. If we can show that this intracel- 
lular network is ergodic, then by observing the long-term 
reaction dynamics in a single cell we can obtain statistical 
information about all the cells at stationarity. 

In this paper we develop a general framework for check- 
ing conditions SCI, SC2 and SC3. However the scope 
of our paper is broader than that. We obtain easily com- 
putable bounds for the statistical moments of the underly- 
ing Markov process and investigate when these moments 
converge with time. We also present conditions for the dis- 
tribution of the process to be light-tailed. The significance 
of these results is mentioned in the next section. 



3 Preliminaries 

In this section we discuss the main results of our paper. In 
particular, we explain how conditions SCI, SC2 and SC3 
can be verified without having to simulate the trajectories 
of the Markov process representing the reaction dynamics. 
Intuitively, these conditions can only hold if the Markov 
process has a low probability of hitting states that have 
a very large size. In our case, the states are vectors in 
M'^ and so we can measure their size by using any norm 
on M"^. The central theme of this paper is to demonstrate 
that for many networks, one can analyze the long-term 
behavior of the reaction dynamics by choosing the right 
norm to measure the state sizes. This right norm has the 



i=l 



Vi\Xi\, 



(5) 



where w is a positive vector in M.'^ satisfying certain con- 
ditions based on the reaction network. We specify these 
conditions in the next section. In the subsequent section 
we show how the vector v can be determined for a large 
class of networks by solving a suitably constructed Linear 
Programming Problem. From now on, we refer to \\-\\v as 
the v-norm on M . 

For any positive integer r, let to^ (i) be the r-th mo- 
ment of ||Xj;j,(i)||„ defined by 



mUt)^K{\\Xx,{t)\\l) 



yes 



(6) 



Using Markov's inequality (see [20]) one can show that 
condition SCI holds if for some r > and some constant 
Cr{xo) we have 



SUpTO^o(i) < Cr{xo). 
t>0 



(7) 



Similarly, condition SC2 holds if for some r > there 
exists a constant C'r such that 



lim sup ni^ 

t— >oo 



(t) < Cr for all xo e 5. 



(8) 



Relations ([7| and ([8]) give uniform and asymptotic upper- 
bounds for the r-th moment of the process (||Xa;o(t)||„)t>o. 
Using these relations we can also obtain uniform and 
asymptotic moment bounds for the Markov process 
{Xxg{t))tyQ representing the reaction dynamics (see 
Corollary 4.3 1. Such moment bound results have appli- 
cations in queuing theory and control theory (see [15] )• In 
Theorem 4.2 we show that for certain values of r, (I?) and 
(Is]) hold and the upper bounds can be easily computed. 

Instead of the r-th moment of the process 
(|jXa;jj(i)||i,)j>o, one can ask if the exponential mo- 
ment of this process is uniformly bounded from above. 
This will happen if for some 7 > we have 



supE(e''ll^="o(')ll 
t>o 



= sup^e^ll^ll"p,„(t,y)<(X). (9) 



t>o _„ 
- yes 



If (|9| holds, then for each t > the distribution Pxgit) 
is light-tailed^ in the sense that its tails are majorized by 
an exponential decay. It shows that all the cumulants of 
the distribution p^g (t) exist, which is an important result 
for the following reason. There is a considerable body of 
research dedicated to estimating the moments of the pro- 
cess {Xx„{t))t>a directly without computing the distribu- 
tion functions Px„ (t). For any integer r > 0, one can easily 
write the differential equations for the dynamics of the first 
r moments. However when the reaction network has non- 
linear interactions, this system of equations is not closed 
for any r. Various moment closure methods (see |16j ) exist 
that specify ways to close these equations artificially by es- 
timating the moments approximately. One such moment 
closure method is the cumulant-neglect method which ig- 
nores the higher order cumulants of the distribution p^o (t) 
for all t > 0. Of course this method is only valid when 



the higher order cumulants exist. This is guaranteed if ^ 
holds. In Proposition 4.4 we give conditions for verifying 

We now come to the question of checking condition SC3 
which says that the process {X.j;g{t))t>o is ergodic. This 
can only happen if the state space S is irreducible, which 
means that all the states are accessible from each other. 
Recall the definition oipx{t,y) from ^. Mathematically, 
we say that S is irreducible if for all x,y € S, we have 
Px{ti,y) > and Py{t2,x) > for some ti,t2 > 0. In 
order to check the irreducibility of S, one has to verify 
that there is no proper subset Si d S, such that once the 
process reaches a state in Si, it stays in Si forever. For re- 
action networks, this can often be easily checked from the 
stoichiometry of the reactions. Assuming irreducibility, 
the ergodicity of the process {Xxg{t))t>o can be estab- 
lished by showing the existence of a norm-lik^ function 
V : S ^ [0, oo) such that the process {V{Xxo{t)))t>o has 
the tendency to decrease whenever it reaches a large value 
(see [26]). In our case, V{x) — \\x\\v plays the role of this 
norm-like function and this allows us to verify ergodicity. 
The precise statement is given in Proposition |4.5| 

Suppose that condition SC3 is satisfied and the process 
{Xxg{t))t>o is ergodic with stationary distribution tt. Us- 
ing Theorem [4.2 1 we can identify functions / : 5 — >■ K for 
which J2yes /(y)'^(y) is finite and 



hm E(/(X 

t—^OO 



.o(i))) = E/(2^Wy) 
yes 



(10) 



holds for any xq £ S (see Proposition 4.6). If / is such a 
function, then the ergodic theorem for Markov processes 
(see [27j) says that 



1 /■* 

lim - / fiX^g{s))ds = V /(y)7r(y) almost surely 

"'0 yes 



(11) 



for any Xq & S. From ( 10 ), we can deduce which moments 



of our underlying Markov process converge to their sta- 



tionary value as time goes to infinity (see Corollary 4.71 



Lastly, we also obtain conditions to check if the stationary 



distribution n is light-tailed (see Proposition 4.8) 



4 General Results 

In this section, we formally present the main results of 
our paper. Let (•, •) be the standard inner product on M''. 
From now on we make the following assumption. 

Assumption 4.1 There is a positive vector v € M'' along 
with positive constants Ci, C2, C3, C4, C5 such that the follow- 
ing conditions hold for all x G S 



K 



^Xk{x){v,C,k) <ci-C2{v,x) (12a) 

fc=i 

K 

and ^\k{x){v,Qkf <C3.+Ci{v,x)+c^{v,x)'^. (12b) 



fe=i 



^A postive function V : cS — >■ R is called norm-like if the set 
{x a S : V(x) < c} is compact for any c > 0. 



Given a reaction network, the procedure to find such a v 
is described in the next section. We now explain how this 
vector V is used. Let {Xxo{t))t>Q be the Markov process 
representing the reaction dynamics. Corresponding to it 
we can define a one- dimensional process {Y{t))t>o as 



Y{t)^\\X,,it)\\,^{v,X,„it)), 



(13) 



where the last equality holds because X^oit) and v are 
positive vectors. In this paper we study the process 
(Ar^o(i))t>o by looking at the process (^(i))t>o, which 
is simpler to study even though it is not a Markov process 
by itself. The dynamics of {Y{t))t>o has two components 
drift and diffusion which have the form X]fc=i ^k{x){v, Cfe) 
and J2k=i ■^kix){v,Ck)^ respectively when X{t) = x. 
Hence ( 12 ) gives upper-bounds for the magnitude of these 
two components. Observe that when the process {Y{t))t>Q 
goes above Ci/c2 then it experiences a negative drift, sug- 
gesting that it will move downwards. This observation is 
crucial for our analysis. 

We are now ready to state the main results of the paper. 
Their proofs are given in the Supplementary Material. In 
all our results, we assume that we are given a reaction net- 
work for which a positive vector v satisfying Assumption 
14. II has been found. 



4.1 Moment bounds 

Our first result establishes that for certain values of r, 
relations ((tI) and ((sl) hold. Let Tmax be the positive number 
given by 



2C2 



if C5 > 
if C5 = 0, 



(14) 



where the constants C2 and C5 are as in Assumption 4.1 



Theorem 4.2 For any positive integer r and xq ^ S let 
rrf^ it) he defined hy ^. If r < r^ax then there exist 
positive constants Cr{xo) and Cr such that ([T]) and ([8]) 
hold. 

The values of the constants Cr{xo) and Cr can be ex- 
plicitly computed using a recursive relationship given in 
the Supplementary Material. For any positive integer 
r, let 'i/'^{xo,t) denote the r-th moment of the process 
{Xxg{t))t>Q at time t. Then 'i^{xo,t) is a tensor of rank 
r whose entry at index ii . . . i,. is given by 

K...^A^O,t) =E{X,,{t)X,,{t) . ..X,^{t)) 



yes 



■ Vir-Pxoit^y) 



for any zi,...,v S {I7 2, . . . ,(i}, where Xxg{t) — 
{Xi{t),...,Xd{t)), y = {yi,...,yd) andpa;o(t) is the dis- 
tribution of Xxg{t). Note that if f = {vi, . . . ,Vd) then 
Xi{t) < \\X.j;o{t)\\v/vi for any i. This gives us the follow- 



ing corollary of Theorem 4.2 



Corollary 4.3 (Moment Bounds) Let 'i'^{xo,t) be the 
r-th moment of the process {Xxg(t))t>o at time t. If 
1^ < ''max then there exist positive constants Cr{xo) and 



Cr such that for any ii, . . . ,ir G {1, 2, . . . , d} 

t>0 [[j=l^h 



have 



and 



O7- 



limsup^^^ j^(a;o,i) < ^fr for all xq e S 



This corollary gives uniform and asymptotic moment 
bounds for the r-th moment of the Markov process 
{Xx^{t))t>Q. Observe that if C5 = then r^ax = 00. In 
this case, Theorem |4.2| says that for each positive integer 
r and xq G 5 there exists a constant Crix^) such that 
([7]) holds. By showing that we have a C > such that 
Cr{xo) < r!C"' for all integers r > 0, we obtain our next 
result. 



Proposition 4.4 Suppose that Assumption 4^.1 holds 
with C5 = 0. Then there exists a 7 > such that 



supE(e''ll^-oWII" 



sup > e 

t>o ~i 
- yes 



invW^ 



Pxo{t,y) < 00. 



This proposition gives sufhcient conditions for checking 
that the distribution p^g (t) is light-tailed for each t > 0. 

4.2 Ergodicity and Moment Convergence 

The next result verifies the ergodicity of a reaction net- 
work satisfying Assumption |4.1[ It follows from Theorem 
7.1 in Meyn and Tweedie [26] . 

Proposition 4.5 (Ergodicity) Let {Xxg{t))t>a be the 
Markov process representing the reaction dynamics and 
assume that the state space S is irreducible. Then this 
process is exponentially ergodic in the sense that there ex- 
ists a unique distribution tt G V^S) along with constants 
B,c > such that for any xq G 5 



sup \p^„ (i, A) 

Acs 



tt{A)\ < Be-"'' for alH > 0. 



This result says that as t — > 00, the distribution Pxo{t) 
converges to tt exponentially fast. Here tt is the unique 
stationary distribution for the dynamics. Henceforth we 
assume that the process {Xxg{t))t>Q is ergodic with sta- 
tionary distribution tt. Using Theorem |4.2| we can prove 
the following. 

Proposition 4.6 Let f : S ^f M. be a function such that 
for some positive integer r < (rmax "1); there exists a 
C > satisfying 

\f{x)\ <C(l + |lx|i;) for all xG 5. 



Then 'Yliy^s f{y)'^{y) is finite and the relations (10 1, (11) 
hold. 

For any positive integer r, let H*" denote the r-th mo- 
ment of the stationary distribution tt. Then H'' is a tensor 
of rank r whose entry at index ii . . . v is given by 

^n...t^ = ^Vii ■■■V^.'^{v) 
yes 



for any ii, . . . , i,. G {1, 2, . . . , d}. Proposition 4.6 allows us 



to prove convergence of moments as time tends to infinity. 



Corollary 4.7 Let '^^'{xo,t) be the r-th moment of the 
process {Xxg{t))t>a at time t. If r < (r„iax — 1) then 
^'■(a;o,i) -^H'' as t^ 00. 



Another consequence of Proposition 4.6 is that for any 
xq G S and positive integer r < (r^ax — 1), the quan- 
tity m]^g{t) (defined by ^) converges as i — > cx) to 
J2yes ll^llXy)- Then Theorem 



4.2 



implies that there ex- 



ists a positive constant Cr such that 



^\\y\\:7r{y)<dr. 



(15) 



yes 



This relation can be used to obtain moment bounds for the 
stationary distribution. Note that if C5 = then r^ax — 00 
and (15) holds for each r. By proving the existence of 



a constant C > such that Cr < r\C^ for all positive 
integers r we get our last result which shows that the 
stationary distribution is light-tailed. 



Proposition 4.8 Suppose that Assumption \41\ holds 
with C5 = 0. Then there exists a 7 > such that 

V"e'''"'^""7r(?;) < 00. 
yes 

The framework described above is very general and can be 
applied to any type of networks that satisfy Assumption 
|4.1[ Below, we specialize the results for two wide classes of 
networks with mass-action kinetics, namely reaction net- 
works with monomolecular and bimolecular reactions. It 
will be, however, pointed out in the the examples that the 
scope is much broader since more general propensities, 
such as those involving Hill functions, can be considered. 



5 Results for affine stochastic re- 
action networks 

Using the analysis tools developed in the previous sections, 
several very general and broad results can be stated for 
the class of affine reaction networks, which is the most 
simple class that can be considered. These networks are, 
in some sense, analogues of linear dynamical systems. Let 
us consider an affine reaction network which involves d 
species that interact through K reaction channels taken 
among the family 





Si 

Si 






Si 


jyj=i'fsj 



(16) 



where i = 1, 



. ,d, A: gN and I'f G Nq. 



The reaction rates 



fcp, k^ and /cf are positive integers. In accordance with 
(H), the reactions are indexed from n = 1 to K, and corre- 
sponding propensities and stoichiometrics are denoted by 
A„(x) and C„, respectively. In the following, it is tacitly 
assumed that the state-space of the underlying Markov 
process describing the stochastic chemical reaction net- 



work ( 16 ) is irreducible. 



5.1 Theoretical results 

Let us start with several theoretical results that character- 
ize the long-term behavior of affine networks of the form 



(16) 



Theorem 5.1 (Nominal ergodicity) Let us consider 
the general affine reaction network ( |16[ ) and assume the 
state-space of the underlying Markov process is irreducible. 

and b e IR>g be further defined 



Let the matrices A G 
as 

K 



fdxd 



^A„(a;)(w,C„> =x''Av + b^i 



(17) 



Then, the following statements are equivalent: 



L The matrix A is Hurwitz, i.e. all its eigenvalues lie 
in the open left half-plane. 

2. There exists a positive vector v G M.'^ such that Av < 
0. 

3. The stochastic reaction network has all its moments 
bounded and converging. 

Moreover, when one of the above statements holds, the 
Markov process describing the reaction network is expo- 
nentially ergodic. o 

The above result shows that, for affine networks, ergod- 
icity and the existence of moment bounds can be directly 
inferred from the properties of the matrix A defined in 



( 17 ). Checking directly whether the matrix A is Hurwitz is 



indeed an easy problem for which efhcient numerical tech- 
niques exist. The second statement, on the other hand, 
characterizes Hurwitzness of A in an implicit way, through 
the existence of a positive vector v > 0. This statement 
may seem superfluous, at first sight, since simply com- 
puting the eigenvalues of A appears to be sufficient for 
our problem. It turns out, however, that it is the most 
important one since finding v > such that Av < is 
a linear programming problem [1]. A formulation that 
is much more flexible than the eigenvalue-based one, es- 
pecially whenever constraints on v needs to be taken into 
account, as this will be the case when analyzing both affine 
and quadratic networks in the rest of the paper. 

Example 5.2 Let us consider the following reaction net- 
work 





^1 



ki 



Si, 

S1 + S2, 



Si 

S2 





S2 + S1 



(18) 



where i — 1,2. Then, the moments go unbounded if and 
only (f7i72 — ^12^21 ^ since the matrix A is not Hurwitz 
in this case. 



We extend now Theorem |5.1| to the more pragmatic case 
of poorly known networks. We therefore assume, in what 
follows, that the structure of network (the reactions and 
stoichiometrics) is exactly known, but that the reaction 
rates are subject to uncertainties. This kind of uncertainty 
reflects naturally into the following parameter-dependent 
matrices: 

n ■n 

A{5) = Ao + ^ 5^A, and b{6) = ba + Y, ^^^^ (19) 



where 5 G [—1, 1]**, rj < K, are normalized time-invariant, 
uncertain and independent parameters. The matrices ^0 
and feo are the nominal matrices which contain the nomi- 
nal values of the rates. The other matrices represent the 
uncertainty radius about the nominal values of the rates. 
We assume further that A{5) and h{5) do not share com- 
mon parameters and that the matrices 

A+ := sup {A{5)} and b+ := sup {b{6)} (20) 



sel-i.i]" 



S^[-l,i\-n 



are well-defined in the componentwise sensq^ This means 
that A_|_ > A{5) for all 5 G [—1,1]'' componentwise, and 
that there exists a i5* G [—1,1]'' for which wc have A^ — 
A(5*). When such matrices exist. Theorem 5.1 admits the 
following robustiflcation: 

Theorem 5.3 (Robust ergodicity) Let us consider 
the general affine reaction network ( |16[ ) described by the 
uncertain matrices ( 19 1 that we assume to the admit 
upper-bounds A^ and b^ defined in (20). Assume further 



that the state-space of the underlying Markov process is 
irreducible for all uncertain parameter values (5 G [—1,1]''. 
Then, the following statements are equivalent: 

1. The matrix A{5) is Hurwitz for all S Cz [—1,1]''. 

2. The matrix j4+ is Hurwitz. 

3. There exists a positive vector w G M'^ such that A^v < 
0. 



4. The stochastic reaction network (16)-(19l has all its 
moments bounded and converging. 

Moreover, when one of the above statements holds, the 
Markov process describing the reaction network is is ro- 
bustly exponentially ergodic. o 

The important point here is that, whenever the matrices 
A^ and 6+ exist, checking ergodicity of a family of net- 
works is not more complicated than checking ergodicity of 
a single network. Note, however, that such a matrix A^ 
may not exist due to the presence of coupled entries in the 
matrix A. This is, for instance, the case when reactions 
of the type Xi — > Xj are involved. In such a case, the 
above result can be conservatively modified to consider 
this particular case. 

5.2 Numerical results 

Several numerical results accompanying the theoretical re- 
sults of the previous section are described in the follow- 
ing. Many computational results can be extracted from 
the previous sections, however, only ergodicity verifica- 
tion and first-order moment bounds computation are ad- 
dressed. The asyrnptptic first-order moment bound, de- 



fined in Theorem 4.2 is given by Ci — ci/c2. So the 
question is, what is the smallest value for such a ratio? 
Or, in other words, what is the smallest compact set that 
eventually contains the first-order moments of {v,X{t))'! 
Several numerical methods will be discussed to solve, ex- 
actly or approximately, this problem. 

The first result is a purely algebraic result which gives 
the optimal (minimal) ratio C]^: 

^Given two real matrices A and B with identical dimensions, then 

J- 



Theorem 5.4 Let us consider the general affine reaction 
network ( 16 ) and assume that the matrix A is irreduciblPl 
Then, we have that 



CI 



h-^VF{A) 
MA) 



where Xf{A) and vp{A) are the Frohenius eigenvalue and 
the Frohenius right-eigenvector of A, respectively. 

When the matrix A is reducible, the Frohenius eigenvec- 
tor vf{A) may contain entries, which would violate the 
conditions for ergodicity that require a positive v. Addi- 
tionally, when V contains entries that are equal to 0, the 
set {x g Ng : {v,x) < C1/C2} is not compact, an irrele- 
vant result when we are interested in finding an attractive 
compact set for the first-order moments. 

The following optimization problem, fully equivalent to 
Theorem \5.1\ can be used in order to deal with reducible 
matrices, by imposing constraints on v. 

Optimization problem 5.5 Let us consider the general 



affine reaction network ( 16 ) and assume that the optimiza- 
tion problem 



s.t. z > Q,v > e 

{zL + A)v<0 



(21) 



is feasible, with [z* ,v*) as minimizer, for some e G IR>o- 
Then, we have 

dl < h^v* jz* (22) 

and Theorem \5.1\ holds. 

Since the feasibility of the above optimization problem 
is equivalent to Theorem |5.1[ this means that when the 
optimization problem is not feasible (or feasible for some 
nonpositive z), then all the moments go unbounded. 

A striking point concerning the above optimization pro- 
gram is that the numbers of variables and constraints are 
given by d -|- 1 and 2d -\- 1, respectively. This means that 
the optimization problem scales linearly with respect to 
the number d of species in the network, and is indepen- 
dent of the number of reactions K. Therefore, from the 
point of view of this optimization problem, the size of an 
affine network can be assimilated to the number of species, 
not the number of reactions. Despite being nonlinear, the 
above optimization problem can be efficiently solved using 
a bisection algorithm over z that is globally and geomet- 
rically converging to z* . Each iteration consists of solving 
a linear program, a class of optimization problems known 
to be very tractable, and for which numerous advanced 
solvers exist [4]. These properties, all together, make the 
overall approach highly scalable, an indispensable prop- 
erty when having the ambition of considering very large 
real- world networks. 



Following (20), the optimization problem (21 1 admits 



the following extension to uncertain networks: 



^A matrix A is irreducible if it is not similar via a permutation 
to a block upper triangular matrix. Note that this has absolutely no 
connection with the irreducibility of the state-space of the Markov 
process describing the network. 



Optimization problem 5.6 Let us consider the general 



( 19 1 which admit the upper-bounds A^ and 6+ defined in 



affine reaction network ( 16 ) with uncertain rate matrices 



(20 1. Assume, further, that the optimization problem 



max^^t, z 

s.t. z > 0,v > e 

{zL + A+)v < 



(23) 



is feasible, with {z*,v*) as minimizer, for some e E lli>o- 
Then, we have 



CiiS) <C; <h\v*/z* 



(24) 



for all 5 E [—1, 1]'' and Theorem 5.3 holds 



Note that when the matrix A_|_ is irreducible, then the 
Perron-Frobenius Theorem can be applied to get an opti- 
mal CI. 



6 Results for quadratic stochastic 
reaction networks 

Similar results are now presented for quadratic stochastic 
reaction networks which, in addition to the affine reactions 



(16), also involve quadratic reactions of the form: 



Si 

Si 



2-^1=1 



ijk , 



(25) 



ijk 



defined for i,j = l,...,d, k g N, v/ G Nq and 
i = l,...,d. The reaction rates fcf^ and A:° are posi- 
tive integers. As for affine networks, it assumed here that 
the state-space of the underlying Markov process is irre- 
ducible. 



6.1 Theoretical results for quadratic net- 
works 



When quadratic reaction networks of the form (16)-(25) 



are considered, the left-hand side of condition (12a) can 
be expressed as 



K 



Xk{x){v, (k) = x''M{v)x -\- x''Av -\- b'^v 



(26) 



i=l 



where M{v) e R'^'"^, A e R'^'"^ and b € R^^. The matrix 
M{v) is a symmetric matrix depending affinely on v. Note 
that, in this case, the matrix A may not be Metzler, al- 
though this will be likely the case for most of the realistic 
reaction networks. 

Let S* := [Ci ... Ck\ be the stoichiometry matrix of 
the quadratic reaction network ( 16 )-( 25 ), and let Sq be the 



restriction of S to quadratic reactions, only.Define further 
the set 

Mq-.^lveR"^ : v> 0, v''Sq = 0} 

which characterizes the set of all positive vectors v such 
that v'' lies in the left null-space of Sq. When v S Mq, the 



quadratic term x'^M(v)x in (26) vanishes, and equality 
( 26 1 reduces to 



K 

E 
1=1 



Xk{x){v,Ck) = x'^Av + b'^v 



which is exactly the same expression as in the case of affine 
networks. This means that, with the additional constraint 
that V g Afq, all the results derived for affine networks 
directly apply to quadratic networks as well. Along these 
lines, we obtain the following result: 

Theorem 6.1 (Quadratic networks) Let us consider 
the quadratic reaction network of the form ( |16[ )-(25) which 
we assume to admit a non-empty J\fq. 

Then, if there exists a vector v G J\fq such that the 
inequality Av < holds, the stochastic quadratic reac- 
tion network (16 1 -(25) is ergodic and has all its moments 
bounded and converging. o 

Unlike for affine networks where v is unconstrained, 
Hurwitzness of A is, henceforth, a necessary condition 
only for Theorem 6.1 to hold. Therefore, the location 



of the eigenvalues of A does not contain all the neces- 
sary information for concluding on ergodicity of the un- 
derlying Markov process governing the reaction network. 
Optimization-based formulations will play a key role in de- 
termining a suitable v that satisfies the conditions above. 

It should be, moreover, pointed out that the condition of 
a non-empty set Afq is central in the statement of Theorem 
|6.1[ When this set is empty, the result cannot be applied in 
any way. A necessary and sufficient condition for J\fq to be 
non-empty is that Sq is rank-deficient, which corresponds 
to a certain form of rarity for the quadratic reactions. This 
condition may seem restrictive at first sight, but it will be 
shown that several important reaction networks from the 
literature actually fall into this category. 

Whenever Afq is empty, linear programming can still be 
applied provided that the matrix M{v) exhibits a partic- 
ular sparsity pattern defined below: 

Definition 6.2 We say that a symmetric real matrix Z = 
[zij] is of class Si if the elements of Z satisfy the condi- 
tions 

for all i ^ j . A network is said to be Si -quadratic if 



its matrix M{v), defined in (26 1, is of class Si for all 

The underlying reason for considering this class of matri- 
ces lies in the fact that, whenever a symmetric real matrix 
Z S R'^^'^ is of class Si, checking whether x'^Zx < for 
all X G M>g is equivalent to the condition that Z < 
componentwise. Based on this fact, the following result is 
obtained: 

Theorem 6.3 (5/-Quadratic networks) Let u: 



sider the quadratic reaction network of the form ( 16 1 



(25) which we assume to be Si-quadratic. Assume further 



that there exists a vector v > such that Av < and 
M{v) < 0. 

Then, the quadratic stochastic reaction network is er- 
godic and has its moments up to order [1 -f 2C2/C5J — 1 
bounded and converging. o 



6.2 Numerical results for quadratic net- 
works 

It is shown here that, once again, the theoretical results 
can be easily turned into linear programs that can be 
checked in a very efficient way. The following result is 
the numerical translation of Theorem 16.11 



Optimization problem 6.4 Let us consider a quadratic 



reaction network (16)-(25) admitting a nonempty set Afq. 



Assume that the optimization problem 

maxz^i, z 

s.t. z > 0,v > e 

{zI + A)v < 
wTS*, = 0. 



(28) 



is feasible with {z*,v*) as minimizer. Then, we have 

CI <b'^v*/z\ V* eAfq (29) 

and Theorem \6.1\ holds. 

Note the presence of the additional constraint v'^ Sq — 
which imposes to v'^ to lie in the left null-space of Sq . We 
can see, again, that the complexity scales linearly with 
the number of species since the number of variables and 
constraints is d -I- 1 and 3d -\- 1, respectively. Therefore, 
quadratic networks admitting a nonempty set Afq may be, 
in fact, not more complicated than affine networks. 

The following optimization problem is the computa- 
tional counterpart of Theorem |6.3| 



Optimization problem 6.5 Let us consider a Si- 
quadratic reaction network of the form (16) -(25). Assume 



further that the nonlinear optimization problem 



max, „ z 



s.t. z > 0,v > e 

{zI-{-A)v < 
M{v) < 0. 

is feasible with {z*,v*) as minimizer. Then, we have 



(30) 



and Theoremw, 



CI < b'^v*/z* 
holds. 



(31) 



The above optimization problem docs not scale as nicely as 



the optimization problem (30) since, in the worst case, the 
number of constraints related to M{v) grows as d{d-\-l)/2, 
hence quadratically in terms of the number of species. 
The number of variables is, however, still equal to d -I- 1. 
Despite the complexity increase, the problem remains 
tractable due to the linear program structure. 

7 Examples 

7.1 Finding an attractive compact set for 
the first-order moments trajectories 

The goal of this section is to describe the computation of 
compact sets that eventually include the first order mo- 
ment of {v,X{t)). Since we are interested in finding the 
smallest one (or an approximation of it), the goal is there- 



fore to minimize the ratio Ci/c2 in (12a), a task which 



can be performed using the optimization problems (21) 



or ( 28 ) , according to the type of networks that is consid- 



ered. Due to the moment closure problem [16 , analytical 
expressions for the steady-state values of the moments of 
quadratic networks are not available. This class of net- 
works is, therefore, the most interesting class to consider. 
Let us indeed consider the following quadratic network 



S^ (32) 




Si + Si 

S2 


k 

fcl2 ^ 


Si, 

S2, 

0. 


Si 

S2 


71 , 
fe2l , 



Si 


72 ^ 






l',{t) 



representing a dimerization process, i.e. Si dimerizes in 
<S'2. It is easily seen that this network is irreducible since 
any point in the state-space can be reached from any 
other point in a finite number of reactions having nonzero 
propensities. Choosing v in A/", where Sq — [—2 l] , 
e.g. v'^ = [l 2] , theoretical calculations show that con- 



dition (12a) holds with c\ — k and C2 — min{7i,72}. This 



proves that the system is exponentially ergodic and that 
all the moments are bounded and converge to a unique 
steady-state value. Solving now the optimization problem 
(28) with numerical values fc = l, 71=72 = 0.2 (the 
values of ki2 and k2i are unimportant here), we get that 
c* = 5, which coincides with the theoretical value. To 
validate the compact set calculation, Monte-Carlo simula- 
tions are performed and yield the values Ei — 1.17 ± 0.01 
and E2 = 1.927 ± 0.02 for the equilibrium values of the 
average number of species Si and 5^2- Evaluating then 
{v, E) with these values yields 



(u, E) =Ei+ 2E2 = 5.024 ± 0.05, 



(33) 



rem 
To 



4.2 



showing that the compact set agrees very well with the 
equilibrium values. Note also that, by virtue of Theo- 
__ we also have limi^oo E[(w,X(t))] < cl / 02- 
illustrate these results, several trajectories for 
E[Xi{t)] and E[X2(t)] for different initial conditions are 
plotted in Fig. 1 where we can see that all the trajectories 
converge to a point, inside the compact set (the surface 
below the dashed line), very close to the boundary. Note, 
moreover, that the trajectories starting from inside the 
compact set never cross the boundary of the set, as pre- 
dicted by the theory. 



7.2 Feedback loop 

Let us consider the feedback loop network of Fig. 2 which 
can be represented by the set of reactions 



Si 

S3 

Q. 


'=2 , e 1 £j 




S2 + S2 


fe23 , 


^1 

^3 


Si 1 ^2, 

k32 , q f q. 


S2 1 •52, 

"'■' s (A 





(34) 
which represents a gene expression network, with Si as 
mRNA and S2 as protein, implementing a feedback loop 
for controlling its own expression through the protein 
dimer S3. The function /(•) can either be an activator 
function (positive feedback) or an inhibitor function (neg- 
ative feedback), but the results are valid for any bounded 
function /(•). Choosing v G Afq and solving for the con- 
ditions of Theorem 6.1 yield that ci — wi sup^>Q{/(a;)} 



Figure 1: Trajectories of the first order moments ^i(i) = 
E[Xi{t)] and ^2(1) = E[X2{t)] of network ^ for differ- 
ent initial conditions (averaging is performed over 5000 
cells). The trajectories converge to the unique steady- 
state value located inside the compact set (the surface 
below the dashed line), very close to the boundary. 



ti 



AAAAA 



/(S3) 



■-> 



■-> 



■-> 



Figure 2: Feedback loop with arbitrary feedback rule. 



and that C2 can be chosen either as C2 = min{72,73} if 
71 > min{72,73} or can be made arbitrarily close to 71, 
otherwise, through a suitable choice for v G Afq. We there- 
fore have the following concluding statement: 

Result 7.1 For any values of the rate parameters and any 
bounded function /(•), the feedback loop with dimerization 
( 34 ) is ergodic and has all its moments bounded and con- 



verging. 

The following stronger result generalizes the above re- 
sult to a feedback loop with multimerization of arbitrary 
order: 

Result 7.2 For any values of the rate parameters and 
any bounded function /(•), the feedback loop with arbitrary 
multimerization and full or partial degradation is ergodic 
and has all its moments bounded and converging. 

7.3 Stochastic switch 

Let us consider the stochastic switch of [33] described by 
the affine stochastic reaction network 







Msi) 



S°i, 



^^^ S°, 



7.., , 



S°i 

S°2 



50 + si 

s° + s^ (35) 



Above S° and Sj represent mRNAs and proteins of gene 
i, respectively. The functions /i(-) and /2(-) are defined as 
inhibiting Hill functions but, as in the example above, any 
bounded functions are allowed. We can state the following 
result: 

Result 7.3 For any values of the rate parameters and any 
bounded functions /i(-) and f2{-), the stochastic switch 



reactions represent how infectious people are recovering 
and how recovered people become susceptible again. 
We then have the following result: 

Result 7.7 For any values of the rate parameters, the 



SIR-model (36) is ergodic and has all its moments hounded 



and converging. 



(35) is ergodic and has all its moments hounded and con- 7.6 Circadian clock 



verging. 

As for the feedback loop, the above result generalizes to 

Result 7.4 For any values of the rate parameters and any 
functions hounded /i(-) and f2{-), the stochastic switch 
with arbitrary multimerization order and full or partial 
degradation is ergodic and has all its moments hounded 
and converging. 

7.4 Repressilator 

We consider here a stochastic repressilator [7j model in- 
volving N genes and negative feedback implemented in 
terms of Hill functions denoted by fi{-), where i corre- 
sponds to the gene number; see Fig. 3. We have the fol- 




Figure 3: iV-gene repressilator. 

lowing result: 

Result 7.5 For any values of the rate parameters and any 
functions fi{-) 's, the stochastic N-gene repressilator is er- 
godic and has all its moments hounded and converging. 

The result can be generalized to account for arbitrary mul- 
timerization of proteins to: 

Result 7.6 For any values of the rate parameters and any 
functions fi{-) 's, the stochastic N-gene repressilator with 
arbitrary multimerization is ergodic and has all its mo- 
ments hounded and converging. 

7.5 Stochastic SIR model 

We consider here the following SIR-model, similar to the 
one in [6], defined as 



Let us consider the circadian oscillator of [33], see Fig. 4, 
which is a network involving 9 species and 18 reactions. 
Using the rate parameters values of [M], we obtain the 
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Figure 4: Circadian clock model of |34) . 

typical oscillatory trajectory depicted in Fig. 5. Applying 
the developed theory on this model, we obtain the follow- 
ing result: 

Result 7.8 For any values of the rate parameters, the 
circadian clock model of 1341 is ergodic and has all its mo- 
ments bounded and converging. 

Using, for instance, the values of [M] and solving for the 
optimization problem (28), we find that ci = 402.5768 
and C2 = 0.1992, and we obtain 



lim E[{v,X{t))] < 2020.778 



(37) 



for the specific v that has been obtained. When averaging 
over a population of 2000 cells, we obtain the sample- 
average trajectories depicted in Fig. 6, where we observe 
convergence to stationary values, even though the sample- 
paths are oscillatory. From ergodicity, we can even state 
that these fixed points for the sample-averages are globally 
attracting. 

7.7 p53 model 

Let us us consider one of the oscillatory p53 models of 
|13) . which is described by the reactions 




/ 
I 



s, 





0, 


R 


R, 


R 



I, 

0, 
s. 





^ 
fc^ 2/ 




^3 


J^ 0, 


Si 

S2 


k, , 


s 


fcs 


T 





0, 
S3, 



Si 
Si 



f (Sl,S3 ) 
/C4 





S1 + S2. 

(38) 



(36) 
where birth and death reactions represent people entering 
and leaving the process, respectively. The only quadratic 
reaction is the contamination reaction which turns one 
susceptible person into an infectious one. The two last 



where Si is the number of p53 molecules, ^2 the number 
of precursor of Mdm2 molecules and Sa the number of 

molecules of Mdm2. The function f(x,y) = — im- 

X -\- kj 

plements a nonlinear feedback on the degradation rate of 
p53. We have the following result: 
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Figure 5: Sample-path of the species of the circadian clock 
model. 
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Figure 6: Tinie evolution of the sample averages of the 
species of the circadian clock model (2000 cells averaging). 



Result 7.9 For any values of the rate parameters, the os- 



cillatory p53 model (381 is ergodic and has all its moments 



hounded and converging. 

7.8 Lotka-Volterra model 

We consider here the stochastic reaction network 







li] 






Si 
Si 



Si 



Si + Ji 





(39) 



which is an open analogue of the deterministic Lotka- 
Volterra system of [H]. The first set of reactions repre- 
sent immigration, the second one reproduction, the third 
one competition due to overpopulation and the last one 
deaths/migrations. We obtain then the following result, 
which is a stochastic analogue of the results in [S] obtained 
in the deterministic setting: 

Theorem 7.10 Let us define T{v) = [i'i7ij] and assume 
that one of the following assumptions hold: 

1. there exists w > such that the matrix V{v) -I- ^{vY 
is positive definite; 

2. there exists w > such that the T{v) + T{vy is copos- 
itiv^ and /3j — 6i < for all i — 1, . . . ,n. 

^A matrix Z is said to be copositive if x^Zx > for all a; > 0. 



and converging. 



C5 



Then, the stochastic reaction network ( 39 ) is ergodic and 
all the moments up to order 1 ^ — 1 are bounded 
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